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A method to compute the scattering solutions of a spinless Salpeter equation (or a Schrodinger 
equation) with a central interaction is presented. This method relies on the 3-dimensional Fourier 
grid Hamiltonian method used to compute bound states. It requires only the evaluation of the 
potential at equally spaced grid points and yields the radial part of the scattering solution at the 
same grid points. It can be easily extended to the case of coupled channel equations and to the case 
of non-local interactions. 



I. INTRODUCTION 

O 

Numerous techniques have been developed to compute the scattering solutions of a Schrodinger equation. Simple 
Runge-Kutta methods can be performed in the case of a local potential, and discretization of the integration domain 
can be used for a non-local interaction Q . All these techniques can be used because the kinetic energy operator can 
be expressed in terms of a derivative operator. This is no longer true in the case of a spinless Salpeter equation for 
)mJ\ \ which the kinetic energy is a complicated square-root operator. 

■ In a previous paper @], we have developed a method to compute the eigenvalues of a spinless Salpeter equation. 
\ This method relies on the fact that the kinetic energy operator is best represented in momentum space, while the 
potential energy is generally given in coordinate space. It requires only the evaluation of the potential at equally 
t-H ' spaced grid points, and yields directly the amplitude of the solution at the same grid points. This method is derived 
from the Fourier grid Hamiltonian method developed to compute the solution of the one-dimensional Schrodinger 
equation, and consequently was called the 3-dimensional Fourier grid Hamiltonian method. It appears very accurate 
and simple to handle. 

' In this paper, we show that 3-dimensional Fourier grid Hamiltonian method can be used to compute the scattering 
y—i \ solutions of a spinless Salpeter equation (or a Schrodinger equation). We focus our attention on the case of a purely 
i central local potential, but the method can also be applied if the potential is non-local, or if couplings exist between 
different channels. Up to our knowledge, this is the first time that the scattering solutions of the spinless Salpeter 
, equation are presented. 

Our method is outlined in Sec. [H] Test applications of the method are presented in Sec. IIIII and a brief summary 
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is given in Sec. 113 
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II. METHOD 
A. Theory 

We assume that the Hamiltonian can be written as the sum of the kinetic energy T and a potential energy operator 
V. The scattering equation is given by 



T + V 



\*) = E\V), (1) 



where T depends only on the square of the relative impulsion p between the particles, V is a local interaction which 
depends on the relative distance, and E is the asymptotic kinetic energy of the two interacting particles. This equation 
is a spinless Salpeter equation if 



f = \Jp 2 + m\ + y p 2 + m\ - mi - m 2 , (2) 

where m\ and mi are the masses of the particles (we use the natural units Ti = c = 1 throughout the text). Equation 
is a Schrodinger equation if 

T=V- with M = mim2 . (3) 
2/j, mi + m.2 
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In configuration space, Eq. is written 

(r\f\r') + (r\V\r') (r'\^)dr"=E{r\^). (4) 



In the following, we only consider the case of a local central potential 

(r\V\r') = V{r)5(r - f') with r = \f |. (5) 
Consequently, the wave function has the following form 

(f|*) = R l (r)Y lm (f) with r = r/r. (6) 
Using the method developed in Ref. |2( , Eq. can be rewritten as 
2 f°° f°° 

-r dr'r'ui(r') dqq 2 T(q 2 ) ji{qr) j^qr 1 ) + V(r)u l (r) = Eu t {r). (7) 



where ui(r) — rRi(r) is the regularized radial function and where functions ji(qr) are spherical Bessel functions. 
Using the following orthogonality relation 

2 f°° 

-xx' I ji(qx)ji(qx')q 2 dq = S(x~x'). (8) 
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One can show that ui{r) oc rji(kr), with k fixed, is a solution of Eq. Q with vanishing potential. The relative energy 
E is then equal to \Jk 2 + m\ + \/k 2 + m 2 — mi — m,2 in the case of a spinless Salpeter equation, and is equal to 
k 2 /(2/j,) in the case of a Schrodinger equation. 

B. Discretization 

In order to compute the scattering solutions of Eq. @ , we replace the continuous variable r by a grid of discrete 
values r% defined by 

r t = iA with i = 0, 1, . . . , N, (9) 

where A is the uniform spacing between the grid points. Regularity at origin ro = imposes ui(ro) = 0. In the 
following, we always consider potential with a finite range linv^oo r V(r) — (the case of scattering by a Coulomb- 
like potential is not considered here). Outside the range of the potential, the solution is a phase shifted free wave 
function. For a value of r/v = NA sufficiently large, we choose to set arbitrarily u;(rjv) = 1 in order to fix the 
normalization of the wave function. 

As explained in Ref. [2| , the spacing A in the configuration space determines the grid spacing Ak in the momentum 
space. Therefore, we have a grid in the configuration space and a corresponding grid in the momentum space 

S7T 

k s = sAk = — with s = 0, 1, . . . , JV. (10) 
If we note V$ = U(r^), the discretization procedure replaces the continuous Eq. J7J) by a matrix equation 

N-l 

[Hij - E Sij] & - -H iN for i = 1, . . . , JV - 1, (11) 



where 



s=l 



N s •> z_< \\NAJ J J \N J Jl \N 



The discrete solution fa of the linear system (|ll|l gives approximately the values of the radial part of the solution of 
Eq. iJTJ at the grid points: fa ~ The phase shift can be computed by using the values of the wave function at 

two points in the region where the potential is vanishing [j| . 
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This method can also be used in the case of a non-local potential and in the case of coupled channels calculations. 
Some details about the implementation of such problems are given in Ref. @. 

Actually, the scattering solution cannot be obtained directly from Eq. JUJ. For instance, in the case of a zero 
angular momentum solution, it is easy to see that Hin — 0. Consequently, we have to set u;(r/v_i) = 1 and to restrict 
the summation in Eq. (fill) to N — 2 (the point ui{r^) cannot be determined). Other normalization problems appear 
for all values of angular momentum. All are due to the discretization procedure as explained below. 

The 3-dimensional Fourier grid Hamiltonian method relies on relation (JSJ. The equivalent discrete orthogonality 
relation on our grid of points is 

o 2 N 

s=l 

One can thus expect that A^'" = §ij for all values of N and I. Actually, the situation is less favorable. In Ref. 0], 
we show that, for I = 0, we have 

Ajf' ,=0) =fy for i,j 1 V !. (14) 



We have verified numerically that 



lim Ajf = % for i,j = l, N-l, (15) 



iv — > 

lim A^' l>1) « Sij for i,j 1 V 1. (16) 

TV— >oo J 

Consequently, the accuracy of this method becomes poorer when I increases; nevertheless for large enough number 
of grid points, very good results can be obtained. 

For scattering problems, it is also interesting to calculate the values of A^^ quantity. One can also expect that 
A'tv '" 1 = <W for all values of N and I. Actually, it is easy to show that 

A^= 0) = 0. (17) 



For other values of I, we have verified numerically that 



lim A^' ¥0) =0 for i = 1, . . . , N and I even, (18) 



iN 

N—>oc 



lim A^' i#0) = 2S lN for i = 1, . . . , N and /odd, (19) 

The simple way to obtain a correct normalization for the solutions, that is to say a value of 1 for the regularized 
radial part of the wave function at the last point of integration, is to solve two different linear systems with respect to 
the parity of the angular momentum. As we shall show in the next section, the following procedures allow to obtain 
accurate solutions of the scattering problem. 

N-2 

^2 [Hij ~ E5 ij] <f>j = -Hi N-i for i = 1, . . . , N - 2 and I even, (20) 
j'=i 

N-l 

[H^ - ESij] <pj = -H iN /2 for i =l, ...,N-1 and I odd, (21) 
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III. NUMERICAL IMPLEMENTATION 
A. Free solutions 

As noted above, solutions of the non-relativistic and semi-relativistic free Eq. (0 can be expressed in terms of 
spherical Bessel functions. It can be easily shown that the vector {ji jo(kri); i = 1, . . . , N — 2}, up to a normalization 
constant, is a solution of the system i|20|) for I = if k = j^s with s = 1, . . . , N — 2. This vector is no longer a solution 
if k is not a integer multiple of . In this case the computed solution matches the solution of the continuous Eq. Q 
everywhere, except near the last point r/v of the domain of integration. This situation is illustrated with Fig.nf° r the 
semi-relativistic free equation. The regularized radial part of the computed solution for a given energy is presented for 
the same value of the spacing, but for different values of rjv- The relative kinetic energy being fixed, different values 
of t* 7v correspond to different values for the parameter s. We can show on this figure that the computed solutions 
differ from the exact solution when s is not an integer. The situation is not modified by a change of the energy. 

If the angular momentum I is different from zero, then the vector {r*j ji{kri)} is not an exact solution of the system 
(12011 or (|2f \ . Nevertheless, a very good approximation of the continuous free solution can be obtained with correct 
values for the parameters A and r^. Again, the computed and the exact solutions can differ strongly near the last 
point rpf. It is worth noting that the differences between the computed and the exact solutions are much less large in 
the non-relativistic case, whatever the value of I. 

In the free case, the phase shift is expected to be zero. Numerically, the phase shift can be determined by using 
two points of the computed solution. If these points are chosen in the region near Tn, the phase shift found can be 
different from zero. On the contrary, when the two points are taken far from rjy, the phase shift value vanishes. 

B. Gaussian potential 

We have tested our method with different finite range potentials in the case of symmetric or asymmetric systems. 
In this section, we shall only present some results obtained with a Gaussian potential 

V(r) = -V e- r2 / a \ (22) 

for two identical particles mi = m,2 — m. 

In the free case, the computed solution can differ strongly from the real solution near the last point r^. This is also 
the case when a potential is turned on. The phase shift can be computed with two values of the numerical solution 
evaluated at two different points r p and r q . If p or q is too close to N, then the value of the phase shift can be very 
bad. Obviously, the two points must be taken in a region where the potential can be neglected with respect to the 
relative kinetic energy. A good procedure to get reliable phase shift is to compute the wave function with a large 
value of tat. Then, the phase shift can be computed with two adjacent points r p and r p _i as a function of the index 
p. By decreasing the value of p from N, the phase shifts will first strongly vary and rapidly reach a stable value, as 
long as r p is large enough to not fall in a region where the potential cannot be neglected. This situation is illustrated 
with Fig. [3 On this figure, the phase shift for two identical semi-relativistic particles with m—1 GeV is plotted as a 
function of r p for two values of the relative energy and for two values of . The Gaussian potential is characterized 
by Vq = 0.5 GeV and a — 10 GeV -1 . It is worth noting that the variation of the phase shifts are much larger for the 
semi-relativistic case than for the non-relativistic case. 

Scattering states have been calculated with our method for two non-relativistic particles interacting with a Gaussian 
potential. In this case, wave functions and phases shifts can also be computed with a great variety of methods. For 
a large range of relative energy and for angular momentum varying from to 4, we have checked that all approaches 
give same results. Within our method, a relative accuracy of at least 10~ 4 for phase shifts can be obtained with 
a grid containing 200-400 points. Obviously, the interest of our method is to compute scattering solutions in the 
semi-relativistic case. 

It is shown in the appendix that a spinless Salpeter equation with a particular separable non-local potential can 
be transformed into a Schrodinger-like non-local equation. In this case, the scattering semi-relativistic equation can 
be solved directly by the Fourier grid Hamiltonian method, or using the equivalent Schrodinger-like form, by usual 
techniques. We have verified, for several values of the parameters and for different values of the relative kinetic energy, 
that all methods give the same results. This yields a direct verification of our approach. 

Finally, we give, in Table [fl the phase shifts for two identical particles interacting via a Gaussian potential as 
a function of the relative kinetic energy. Results have been computed for a non-relativistic and a semi-relativistic 
kinematics, and for two values of the angular momentum I. As expected, phase shifts are similar for low relative 
energy, and differ when energy increases. 



5 



IV. SUMMARY 

The 3-dimensional Fourier grid Hamiltonian method, used in a previous work to compute bound states 0, appears 
as a convenient method to find the scattering solutions of a spinless Salpeter equation (or a Schrodinger equation). 
It has the advantage of simplicity since it requires only the evaluation of the potential at some grid points and it 
generates directly the values of the radial part of the wave function at the same grid points. Moreover, the method 
can be extended to the cases of non-local interaction or coupled channel equations. Up to our knowledge, this is the 
first time that the scattering solutions of the spinless Salpeter equation are presented. 

Meson-meson scattering has been recently investigated in terms of quark degrees of freedom within the framework 
of the non-relativistic resonating group method From this work, it appears that the use of a semi-relativistic 
kinematics is necessary to avoid inconsistencies related to the non-relativistic formalism. We have calculated a semi- 
relativistic version of the pion-pion scattering equation. This equation is a scattering spinless Salpeter equation. In 
this framework, the method presented here appears suitable to calculate the corresponding phase shifts [3 ■ 

The accuracy of the solutions of the numerical method presented here can easily be controlled since it depends only 
on two parameters: The number of grid points and the largest value of the radial distance considered to perform the 
calculation. This distance must be large enough to fall in the region where the potential can be neglected with respect 
to the asymptotic kinetic energy. Both parameters can be automatically increased until a convergence is reached for 
phase shifts. 

The method involves the use of matrices of order (N x N), where N is the number of grid points. Generally, 
the most time consuming part of the method is the solution of the linear system. This is not a problem for modern 
computers, even for PC stations. Moreover, several powerful techniques exist and can be used at the best convenience. 
A demonstration program is available via anonymous FTP on umhsp02.umh.ac.be/pub/ftp_pnt/. 
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PARTICULAR CASE OF THE SPINLESS SALPETER EQUATION 

The spinless Salpeter equation for two identical particles interacting via a non-local potential can be written 



2v/p 2 + m 2 *(r) =M*(f)- / df'W{r,r')^(r'). (23) 

Acting on both sides with the square-root operator gives 

4(p 2 + m 2 ) *(r) = M [hl^{r) - J dr*' W{r,f')9(r') \ -2 J dr' ^/p 2 + m 2 W{r, r') *(f')- (24) 

If the potential has the form W{r, r') — Vq V(r) V(r'), and if we perform the integrations on angular variables, we 
obtain 

(p 2 + m 2 - R (r) = -irVa (m + 2 v / p 2 + m 2 ) V(r) J™ dr' r' 2 V(r') R (r'). (25) 

where Ro{r) is the radial part of the S-wave function W(r ). It has been shown in Ref. y| that 

2 + m 2 e~ mr = — K (mr), (26) 

where Kq{x) is a modified Bessel function 0, p. 952]. In this case, if we choose V(r) = exp(— mr), then Eq. ill'MI) 
reduces to a non-local Schrodinger-like equation 

d 2 M 2 \ r°° 

— - to 2 + — J «o(r) = V (Mire- mr + 8mK (mr)) r J dr' r' e~ mr u (r'), (27) 
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TABLE I: Phase shifts for two identical particles with m = 1 GeV as a function of the relative kinetic energy E. The interaction 
is a Gaussian potential with Vo = 0.1 GeV and a = 5 GeV -1 (see Eq. 1221 '). Results are given for a non-relativistic (NR) and 
a semi-relativistic (SR) kinematics, and for two values of the angular momentum I. 

I E (GeV) £nr (rad) Ssr (rad) 



0.001 


-0.192 


-0.189 


0.01 


-0.627 


-0.619 


0.1 


1.363 


1.376 


1. 


0.447 


0.524 


10. 


0.156 


0.256 



0.001 


-0.231 


-0.226 


0.01 


-0.931 


-0.925 


0.1 


1.241 


1.254 


1. 


0.479 


0.517 


10. 


0.156 


0.256 



where Uo(r') = r' Ro(r'). 
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FIG. 1: Regularized radial part Ro(r) of computed solutions (dotted lines) and the exact solution (solid line) for the semi- 
relativistic free equation with a given relative kinetic energy and a zero angular momentum. The computed solutions are given 
for the same spacing, but for different values of rjv. The value of the parameter s = NAk/ir is also presented. The energy is 
chosen in order that s — 3 for rjv = 40 GeV -1 . All computed wave functions are normalized to match the exact solution. 
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FIG. 2: Phase shifts 8 in radians for two identical semi-relativistic particles with m = 1 GeV, interacting via a Gaussian 
potential with Vo = 0.5 GeV and a = 10 GeV -1 . The phase shifts are plotted as a function of r p (see Sec. lIIIB^ for two values 
of the relative energy and for two values of rjv. The potential as a function of r is also indicated. 



